Typical and large-deviation properties of minimum-energy paths on disordered 

hierarchical lattices 



0. MelcherlQ and A. K. HartmanrQ 

Institut fur Physik, Universitat Oldenburg, 
Carl-von-Ossietzky Strasse, 26111 Oldenburg, Germany 

(Dated: February 27, 2013) 

We perform numerical simulations to study the optimal path problem on disordered hierarchi- 
cal graphs with effective dimension d e a ~ 2.32. Therein, edge energies are drawn from a disorder 
distribution that allows for positive and negative energies. This induces a behavior which is funda- 
mentally different from the case where all energies are positive, only. Upon changing the subtleties of 
the distribution, the scaling of the minimum energy path length exhibits a transition from self-affinc 
to self-similar. We analyze the precise scaling of the path length and the associated ground-state 
energy fluctuations in the vincinity of the disorder critical point, using a decimation procedure for 
huge graphs. Further, using an importance sampling procedure in the disorder we compute the 
negative-energy tails of the ground-state energy distribution up to 12 standard deviations away 
from its mean. We find that the asymptotic behavior of the negative-energy tail is in agreement 
with a Tracy- Widom distribution. Further, the characteristic scaling of the tail can be related to 
the ground-state energy flucutations, similar as for the directed polymer in a random medium. 
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I. INTRODUCTION 



Many problems in physics and computer science can 
conveniently be modeled using graphs. Thereby it is of- 
ten inevitable to assign attributes to the edges that as- 
sist in specifying the problem under consideration. E.g., 
weighted graphs, where a weight, is associated with each 
edge, might be used to model disordered environments. 
For a given weighted graph the minimum- energy path 
(MWP) problem refers to the paradigmatic optimization 
problem of finding a simple (i.e. loopless) path, connect- 
ing two distinguished nodes of the graph, along which 
the sum of the edge weights is minimal. The MWP prob- 
lem quite naturally lends itself to study a multitude of 
lattice-path models in the context of disordered systems. 
In this regard, it has proven to be useful in order to 
characterize, e.g., linear polymers in random media 0- 
5], domain wall excitations in disordered environments 
such as spin glasses [fUl] and the solid-on-solid model 
@. Due to this relation to physical problems, the weight 
will be denoted as energy in the following. If the disorder 
is drawn from a distribution that allows for nonegative 
edge energy only, as for the canonical "directed polymer 
in a random medium" (DPRM), the groundstate config- 
uration of the polymer can be computed efficiently using 
Dijkstra's algorithm [l(| [HI- However, if the disorder 
distribution allows for edge-energies of either sign, as for 
the problem of finding a minimum energy domain wall 
in 2D Ising spin glasses @, H, [H| (given that there is no 



'Electronic address: 
t Electronic address: 



Oliver . melchert @ uni-oldenbur g . de 



alexander.hartmanniauni-oldenburg.de 



closed path with a negative energy) or more generally for 
the negative-weight percolation (NWP) problem [l3l - fl8l j . 
the solution of the MWP problem requires a nontrivial 
transformation to an auxiliary minimum-weight perfect 
matching problem 19]. Furthermore, the properties of 
MWPs with negative edges are fundamentally different 
from the case where all edge energies are non-negative 

mm- 

To specify the NWP more precisely, one considers, say, 
a regular d = 2 square lattice graph with side length L 
and free boundaries in one direction, periodic boundaries 
in the remaining direction, and energies drawn from a 
distribution that allows for edge energies of either sign. 
The details of the energy distribution are controlled by 
a tunable disorder parameter. For a given realization of 
the disorder one might be interested in, say, an agent 
"harvesting" the negative energies (seeing as a negative 
cost, i.e., a resource, e.g., an energy) along a freely ad- 
justable path between two given points. This means 
the walker might have to use edges with positives en- 
ergies as well, i.e., spend some amount of the resource. 
In addition some resources might be harvested, in par- 
allel or in competition to the walker, by other walkers 
which are not restricted to walk between two given end- 
points. These other walkers are only present for walks 
where the amount of the harvested resource is larger 
than the amount of the resource spent. Mathematically 
this means we consider configurations consisting of a sin- 
gle path and a set of loops, i.e. closed paths, such that 
the total sum of the energies assigned to the edges that 
build up the path and the loops attains a minimum. As 
an additional optimization constraint the path might be 
forced to span the lattice along the direction with the 
free boundaries. This means, the walker covers a large 



2 




FIG. 1: Samples of minimum-energy configurations consist- 
ing of one path (forced to span the lattice along the direction 
with the free boundaries) and a set of loops for a 2D square 
lattice with side length L = 64 and periodic boundaries in the 
horizontal direction. The snapshots relate to different values 
of the disorder parameter p, where (a) p < p c , (b) p ~ p c , 
and, (c) p > p c . In the limit of large system sizes and above 
the critical point p c , paths might span the lattice along the 
direction with the periodic boundaries. 



fraction of the lattice which allows him to maximize the 
amount of the harvested resource. Further the path and 
the loops are simple and are not allowed to intersect each 
other. Therefore, they exhibit an "excluded volume" 
quite similar to usual self avoiding walks (SAWs) [20| . 
A pivotal observation is that the NWP model features a 
disorder driven, geometric phase transition, signaled by 
the emergence of paths that span the lattice along the 
direction with the periodic boundary conditions. In this 
regard, depending on the disorder parameter, one can 
identify two distinct scaling regimes: (i) a phase where 
the paths tend to be short in length, displaying a self- 
affine scaling with system size, see Fig. QTa), and, (ii) 
a phase where the paths tend to be long and exhibit a 
self-similar scaling, see Figs. [TJb),(c). From the previous 
analyses for d = 2 we found that right at the critical 
point, the paths are self-similar with a fractal dimension 
d f = 1.268(1), see Refs. [HEl- 

Here, we study a particular MWP problem (which 
closely resembles the NWP problem) in a Migdal- 
Kadanoff-like renormalization group scheme on hierar- 
chical lattice graphs (constructed using a "Wheatstone 
bridge" elementary c ell, see Fig. [5]) with an effective di- 
mension <i c ff w 2.32 [21H23l |. where an exact decimation 
procedure can be used to analyze huge graphs [24|, |25| ■ 
We address the critical behavior of the MWP in this setup 
by monitoring observables related to the path energy and 
path length. The subtleties of the construction procedure 
that leads to hierarchical lattices with effective dimen- 
sion <i c ff 2.32 even allows to probe the transition from 
the self-affine to self-similar scaling of the path length, 
as observed for the NWP problem on hypercubic lattice 
graphs [HE!. 

Similar to previous studies of minimum- energy path 
problems on hierarchical grap hs [H , [(| [5l|, [26Tj3l| and reg- 
ular lattices [l|, @, la Hfl, l22l 123, we here consider the 
finite-size scaling of the length and energy of the paths 
as well as the associated energy fluctuations. We fur- 
ther complement the simple sampling (SiSa) estimates 
of the path-energy distribution by an importance sam- 
pling (ImSa) procedure in the disorder [3Jt38|, allowing 



to resolve the respective distributions up to 12 standard 
deviations away from its mean. In this regard it is found 
that the asymptotic behavior of the negative-energy tail 
is in agreement with a Tracy- Widom distribution. Fur- 
ther, the characteristic scaling of the tail can be related 
to the path energy fluctuations, similar as for the di- 
rected polymer in a random medium (27j . Note that, 
apart from the disorder distribution, the MWP problem 
considered here is similar to the optimal path problem on 
hierarchical lattices as studied in Ref. [21| (therein, the 
authors considered a uniform distribution of nonegative 
edge-energies, only). 

Above we pointed out that the MWP problem studied 
here closely resembles the NWP problem. At this point 
we would like to point out the major similarities and 
differences of both models: similar to the NWP problem, 
in the MWP problem the sum of energies of edges that 
build up a path is object to minimization. As a major 
difference note that the path found in the context of the 
NWP problem is not necessarily the (absolute) minimum- 
energy path. The reason is that in the NWP problem, 
a global minimum of the energy for a single path plus a 
(possibly empty) set of loops (all with negative energy) 
is searched for, see Fig. [1] This is in contrast to the 
MWP problem on hierarchical lattices, where a particular 
decimation scheme, see sect. HU allows to obtain a truly 
minimum-energy path in a framework where no loops are 
considered. However, also note that at the critical point 
of the NWP model, see Fig.QTb), the "additional" loops 
are small and resemble a rather dilute "gas" of loops 
which arc unlikely to affect the statistics of the path that 
spans the system in between the free boundaries. Hence, 
in the vicinity of the critical point of the model we expect 
the MWP problem studied here to provide a reasonable 
approximation to the path in the NWP problem. 

The remainder of the presented article is organized as 
follows. In section UU we explain the construction pro- 
cedure to obtain the hierarchical lattice graphs and we 
outline the pool method used compute the properties of 
the paths for huge graphs. In section Hill we list the re- 
sults of our numerical simulations in terms of which we 
locate the self-affine to self-similar transition of the path 
length and where we put under scrutiny the path-energy 
distribution. In section HVl we conclude with a summary. 



II. MODEL AND ALGORITHM 

The hierarchical lattices considered in the remainder 
of the presented article can be constructed using a 
simple deterministic rule. This rule specifies how 
the individual edges in a graph at a given iteration 
step / need to be transformed in order to obtain a 
graph at iteration step I + 1. Let Gj = (V/,£r) 
denote a hierarchical graph that consists of a set of 
nodes i e V; and a set £/ C vj 2 ' of undirected edges 
e = {i,j} G £/. The number of edges is given by 
Mi = |£f |- The transformation in order to proceed from 
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FIG. 2: Illustration of the first two iteration steps in the 
construction procedure to obtain the hierarchical graphs for 
which the presented study is carried out. The linear extend 
Li of the graphs Gi is also indicated. As explained in the text, 
the resulting graphs have an effective dimension d e s « 2.32. 



decimation 




FIG. 3: Illustration of the decimation procedure used to fill 
the pools. Repeatedly, five edges ei . . . es are randomly picked 
from pool Vi in order to construct a five-edge subgraph G' 
which is in turn decimated to a single edge e (see text) and 
added to pool Pi+i. 



Gi to Gj+i reads as follows: each edge e = {i,j} G £i 
is replaced by a subgraph G' consisting of four nodes 
{i, k^\ k^\j} (therefore the set of nodes needs to be 
amended by two nodes fcg and fe{ ) and five edges 

{{*, {i, k[% {k^\ {k^,j}, {k[ ij \j}}. 

The nodes i and j are referred to as the terminal nodes 
of the subgraph. After the transformation is completed 
the number of edges increased to M/ +1 = 5 x Mj, 
and the linear extension of the graph has doubled, i.e. 
Lj+i = 2 x Lj. At Z = the construction procedure 
is started with a single edge, meaning that M = 1 
and Lq = 1. Hence, Mj = and L/ = 2 7 . Prom 
the increase of the number of edges Mj as a function 
of the linear extension Lj of the graphs according to 
Mi = 2 /lo S2(5) = L d f s it is possible to obtain the 
effective (fractal) dimension of the hierarchical lattices 
as d e g — log 2 (5) ~ 2.32. The construction procedure is 
illustrated in Fig. [31 where, starting with a single edge 
at I = 0, the two steps Go — > G\ —> G2 are shown 
explicitly. Finally, a path is represented by an ordered 
set of edges. E.g., regarding the subgraph G', a possible 
path that connects its terminal nodes i and j reads 
p=({i,k^h{ki ij \k^},{4 ij \j}). 

The minimum-energy path problem we address here 
reads as follows. Let s, t denote the endnodes of the 
single edge at / = 0. Perform a number of / max itera- 
tion steps to yield a hierarchical graph G/ max and assign 
a random energy to each edge, drawn from a given dis- 
order distribution. Finally, compute a minimum energy 
s-t path for the graph G/ max . The precise topology of 
the resulting path depends on the particular realization 
of the disorder and has length £ G [2 /max , 3 /max ]. Bear in 
mind that in order to compute one such path, a graph 
with M/ max = 5 /max edges needs to be constructed. Con- 
sequently, a number of M/ max random deviates need to 
be drawn from the disorder distribution. A more efficient 
way to sample minimum energy s-t paths for the case of 
hierarchical graphs at large values of I is provided by 
the pool method [24, 25]. Therein one maintains a set of 
/ max pools Vj of (effective) edges, with 1 = 0... 7 max — 1. 
The number of edges in each pool is the same and is de- 
noted by N. An individual edge carries two attributes 



e = (E, £), where E denotes the energy and £ the length 
of a path associated with the edge. At I = the edges 
are initialized with £ = 1 and the energies E are drawn 
from a specified disorder distribution Pq(E), signifying 
the "single edge level" . In order to proceed from pool Vi 
to 7-V+i, the following three-step decimation procedure, 
sketched in Fig. [31 has to be repeated N times: 

(i) pick five edges e\ . . .e§ at random from pool Vi. 
Combine these to form a subgraph G' as explained 
earlier. 

(ii) from the four distinct paths that connect the 
terminal nodes of G', determine the mini- 
mal energy path p*, i.e. the path p* 6 
{(ei,e 2 ), (e 3 ,e 4 ), (ei,e 5 ,e 4 ), (e 3 ,e 5 ,e 2 )} for which 
E* = J2e£ P * E(e) = min. Correspondingly, the 
length of the path reads £* = X)ee P * ^"( e )- 

(hi) set up a new edge having attributes e = (E*,£*) 
and add it to pool Vi+i- 

After Vq has been initialized, all pools up to I = I max 
might be filled in this manner. Note that an edge e G Vi 
effectively corresponds to a hierarchical graph G/, i.e. 
it has a "hidden" substructure that allows to represent 
a minimum-energy path with length t G [2 / ,3 7 ]. The 
attributes of the edge encode the characteristics of the 
respective path, i.e. its energy E and length £. A pool 
Vi thus consists of N instances of minimum-energy paths 
for hierarchical graphs at iteration step /. Thereby, the 
computational resources needed to fill a pool stay con- 
stant as I increases. 

Further, each pool specifies its own distributions Pi(E) 
and Pi{£) of path energies and path lengths, respectively. 
These allow to quantify the scaling behavior of the aver- 
age path length with system size as {£) cx Lj f , defining 
the fractal dimension df of the paths, the average path 
energy (E) cx , and the fluctuation of the path en- 
ergies as var(.E) = (E 2 ) - (E) 2 - L] n . Note that these 
energy fluctuations are measured with respect to the lin- 
ear extension of the considered lattice graphs. This can 
be compared the ground-state energy fluctuations of the 
DPRM, which are commonly measured for polymers of a 



4 



0.94 




0.84 - 



6 8 10 12 14 16 18 20 22 
I 

FIG. 4: Results for the estimation of the critical point p c 
using the secant method. The main plot shows the scaling of 
the iteration step dependent critical points p c (I)- The dashed 
line is a fit to the scaling form p c = p c + a2~ Ib in the interval 
I 6 [5,22], resulting in the estimates p c = 0.8436(2), a — 
O(l), and b = 0.866(4). The inset illustrates the probability 
density function P(p c a) of the effective critical values p e ff for 
/ = 6, 8, 22 where a number of 256 independent estimates 
were considered. 

fixed length L, giving rise to the fluctuation exponent u> 
defined as varies) cx L 2u , see Ref. Q. Hence, to com- 
pare our results with the DPRM case, one should rewrite 
the energy fluctuations as a function of the average path 
length (£). This leads to a corresponding estimate of w 
via using the relation u> — Q/df. 

In the following section, we will use the algorithmic 
procedure outlined above in order to study the minimum- 
energy path problem for an increasing number of itera- 
tion steps / and for a large range of values of the disorder 
parameter p. 



III. RESULTS 

In the presented study, the disorder distribution is a 
Gaussian with mean p and width a 2 = 1. A tunable 
disorder parameter is defined as p = 1/p, so that the 
standard normal distribution is recovered in the limit 
p — > oo. Typical values for the pool-size and iteration 
steps are N = 10 6 and / max = 30, respectively. 

To facilitate intuition, note that as pi — > oo, typ- 
ical minimum-energy paths will exhibit positive ener- 
gies. Thus, an increasing path length will lead to an 
increasing path energy. Hence, a minimum-energy path 
will tend to be short in length. Considering hierarchi- 
cal lattice graphs Gi, one might consequently expect 
a scaling behavior (£) oc 2 1 , implying a scaling expo- 
nent df = 1. On the other hand, as p — > — oo, typi- 
cal minimum-energy paths will exhibit negative energies. 
Therefore, an increasing path length results in a decreas- 
ing path energy, leading to expect {£) oc 3 7 , and therefore 
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FIG. 5: Variability of the energy and length of the minimum- 
energy paths contained in pools corresponding to different 
iteration steps /. The data is normalized in a way that curves 
for different values of / overlap as p — > 0. (a) The main plot 
shows the variance of the path energies, where the curves 
indicate a change in the scaling behavior at p m 0.8367. This 
is supported by the scaling of the average path energy, shown 
in the inset, (b) The main plot shows the variability of the 
path length, and the inset illustrates the scaling behavior of 
the effective critical points p e g(I) that indicate the associated 
peak position. 



df = log 2 (3) w 1.585. In between these two extremal 
"trivial" cases, there exists a particular value p c = 1/ p c 
of the disorder parameter that signifies the onset of "pro- 
liferation", where (E) — as / — > oo. Here the average 
path length exhibits a non-trivial scaling behavior dis- 
played by a scaling exponent 1 < df < 1.585. 

At first we attempt to estimate the value of p c by 
means of the secant method [3^], considering different 
initial pools at a given value of /. Next, we consider one 
particular initial pool to quantify the scaling behavior 
of the average minimum-energy path length and energy. 
Finally, we put under scrutiny the probability density of 
minimum-energy path energies. 
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FIG. 6: Critical exponents that characterize the minimum- 
energy path length and energy on hierarchical lattice graphs 
for different values of p below, right at, and above the effective 
critical point p e g — 0.836688(f) (computed for one exemplary 
pool of size N = I0 6 considering / < 30). The main plot 
shows the fractal dimension df obtained from a power law fit 
to the scaling form (£) oc 2 Id f for / £ [20, 30] (only exception: 
the fit at the critical point p c was restricted to the interval 
I 6 [5 : 10]). The inset shows the energy fluctuation exponent 
fi obtained from a fit to vax(E) oc 2 2m in a similar manner. 



A. Location of the critical point where the 
minimum-energy path energy vanishes 

In order to approximate the critical point p c we con- 
sidered pools of size N = 10 5 and I < 22. So as to 
arrive at an estimate of p c at iteration step I we pro- 
ceeded as follows: Using the secant method we prepared 
a number of M = 32 independent estimates of effec- 
tive critical values p^g, where (E) w at the considered 
value of /. The distribution of these effective critical 
values (see inset of Fig. 2]) is characterized by the aver- 
age p c (I) = (1/M)££ipS- E.g., at I = 22 we yield 
p c (I) = 0.8435(7), wherein the standard deviation among 
the M independent estimates reads <j Pc = 0.004. The av- 
erages exhibit the scaling behavior p c (I) = p c + a2~ Ib , 
where a fit [4(| to the interval / £ [5, 22] yields p c = 
0.8436(2), a = 0(1), and b = 0.866(4), see Fig. H 
The results did not depend much on the pool size, e.g. 
considering N = I0 4 and proceeding as above we find 
p c = 0.8435(1) and b = 0.868(5). 



B. Trivial to non-trivial transition of the average 
minimum-energy path length 

As it appears, for large values of I and p c f» 0.84 one 
should observe a vanishing average path energy. In the 
presented subsection we consider a single pool of size 
N = 10 6 and I < 30 in order to assess the scaling behav- 
ior of the minimum-energy path length and energy. For 



that pool we find that the average path energy changes 
its sign at p « 0.8367 (see inset of Fig.JSJa)). Further, at 
that approximate value the scaling behavior of the fluc- 
tuations related to the path energy and length change 
significantly (see Figs. [5ja),(b)). In this regard, the pre- 
cise location of the peak position related to var(£) (see 
Fig. [HJb)) can be used to define an iteration-step depen- 
dent effective critical point p c s(I). Seen as a function 
of /, these effective critical points can be used to pin- 
point the precise location where the proliferation tran- 
sition of the path length occurs in the limit I — > oo. 
The effective critical values exhibit a scaling of the form 
Pcs(I) = p c +o,2~ Ib ', where a fit to the interval I £ [f 0, 20] 
yields the estimates p c = 0.836688(f), a — 0(1), and 
b = 0.76(2), see inset of Fig. [5fb) . As pointed out above, 
for p < p c an increase in path length most likely results 
in an increasing path energy. Hence, one can expect that 
for p < p c the minimum-energy path problem investi- 
gated here effectively corresponds to the optimal path 
problem studied in Ref. [21], wherein a nonegative uni- 
form disorder distribution was considered. From this it 
is immediate to expect df = 1, cIe — 1 arid Q = 0.3 for 
p < p c - From a direct fit to the scaling forms (£) oc 2 Idf , 
(E) oc 2 IdE and var(E) oc 2 2ln we obtained the numer- 
ical values for the scaling exponents df, ds and Q (and 
consequently the "corrected" exponent ui) as shown in 
Fig. [5] and listed in Tab. [T] For values of p below and 
above the critical point, the fits were restricted to large 
iteration steps, i.e. / 6 [20,30]. However, note that it is 
difficult to prepare a system right at p c : as I increases, 
fluctuations will eventually cause the system to assume 
the asymptotic scaling behavior characteristic for p < p c 
or p > p c . Hence, in order to obtain the scaling expo- 
nents for p p c the fitting procedure was restricted to 
intermediate iteration steps I € [5, 10] (I € [5, 15] in case 
of (E)), only. Anyway, for the case of the energy, we 
actually require (E) ~ at p — p c , hence the value of 
dE right at the critical point is of limited relevance and 
a pure numerical artifact. 



TABLE I: Critical exponents that characterize the self-affine 
to self-similar transition of the minimum-energy path length 
on hierarchical lattice graphs. The table lists the numerical 
values of the critical exponents below, right at, and above 
the effective critical point p c g — 0.836688(1), computed for 
one examplary pool of size N = 10 6 considering I < 30. 
From left to right: numerical values of the scaling exponents 
df and dE for the path length and energy, respectively, the 
energy fluctuation exponent f2 (where the energy fluctuations 
are considered as a function of the linear extension of the 
lattice graph), and the corrected energy fluctuation exponent 
u) (where the energy fluctuations are considered relative to 
the actual path length). 
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FIG. 7: Local scaling exponents obtained for the length and 
the energy fluctuations of the minimum-energy paths, close to 
the critical point p c — 0.836688(1). The main plot shows the 
behavior of the local exponents d\ oc (p) associated to the path 
length for an increasing number of decimation steps / as a 
function of the disorder parameter p, and the inset illustrates 
the scaling exponents f2ioc(p) related to the fluctuation of the 
path energies. 

Another means to quantify the scaling behavior of 
the above observables is given by the local scaling ex- 
ponents. E.g., denoting the average minimum- energy 
path length at a given value of p and iteration step I 
as (£(p))i, one can obtain the local analog to the fractal 

dimension as d[2(p) = log 2 ((£(p))i+i/{i(p))i)- The re- 
spective error can be obtained via error propagation as 

Hoc = mp))i/{Kp))l + S(i(p)) I+1 /(l(p))i+i. The re- 
suiting local scaling exponents are shown in Fig. [7] The 

local equivalent fij^ of the energy fluctuation exponent 
can be computed in similar manner, see inset of Fig. [7] 
For a given value of p and for increasing /, the local ex- 
ponents become independent of / and tend to a limiting 
value that is within error bars in agreement with the nu- 
merical values of df and listed in Tab. |TJ 

C. Importance sampling results for the 
ground-state energy distribution 

As stressed above, for a disorder parameter p < p c 
the minimum energy path problem considered here effec- 
tively corresponds to the optimal path problem studied 
in Ref . [2l| . This is further highlighted by the probability 
distribution function (pdf) of minimum-energy path en- 
ergies. In this regard, Fig. [5] shows the simple-sampling 
estimate of the ground-state energy distribution for the 
minimum-energy path at p = 0.830 for different iteration 
steps /, obtained from pools of size N = 10 6 . As evi- 
dent from Fig. [5J data curves corresponding to different 
iteration steps / scale according to 

P I (E)=a- 1 P{{E-{E))/a E ), (1) 
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FIG. 8: Simple-sampling estimate of the ground-state energy 
distribution for the minimum-energy path at p — 0.830. The 
distribution of the optimal path energy for the case where 
edge-energies are drawn uniformly from the interval [0, 1] (i.e. 
the case considered in Ref. [2l]]) is also shown (in the key, 
the respective symbol is marked as (u)). The dashed line 
indicates the guiding function fitted to the negative tail of 
the ground-state energy distribution. 



where e = (E — (E))/oe defines a reduced energy with 
(E) and o\e describing the average and standard devia- 
tion of the distribution Pj(E), respectively. This means 
that the scaling function P(e) does not depend on the 
value of /. Note that the rescaled distribution of the 
optimal path energy for the case were edge-energies are 
drawn uniformly from the interval [0, 1] (i.e. the case con- 
sidered in Ref. [21| ) assumes the same scaling form and 
falls onto the same master curve. 

Due to the close correspondence between the 
minimum-energy path problem studied here and DPRM 
problem one might expect that the scaling function 
P(e) has the shape of a Tracy- Widom distribution (see 
Refs. [H The ne § ative tail of tne Tracy- Widom 

distribution exhibits the asymptotic scaling Ptw(^) °< 
cxp(— c|a;|' ? ), where in case of the d = 1 DPRM (having 
one space and one time direction) one has t^^ m = 3/2 
exactly [37j . The exponent r), describing the negative tail 
of the pdf P(e), is thereby related to the energy fluctu- 
ation exponent u> (listed in Tab. [TJ) by means of the ex- 
pression rj = 1/(1 — uj). Furthermore, SAWs in quenched 
random hierarchical environments at the critical point 
were studied [27| using real-space renormalization tech- 
niques (quite similar to the approach presented here). 
Also in that case both tails of the scaling function P(e) 
are consistent with an exponential decay as given above. 
Under the assumption that the tails of the distribution 
reproduce under rescaling, they arrive at the estimates 
r/_ = 1/(1 — Q/df) and 77+ = 1/(1 — 17) for the negative 
and positive tail, respectively. However, their numeri- 
cal results then did not allow to conclude with precise 
numerical estimates for the exponents rj±. 

Subsequently we address the question whether the dis- 
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tribution of ground-state energies in the minimum- energy 
path problem for p < p c , p s» p c , and p > p c is consistent 
with a Tracy- Widom scaling form and we attempt to ob- 
tain a numerically precise estimate of the negative tail 
exponent r] for the above three cases. To this end, we 
consider an importance sampling procedure in the dis- 
order [34|, [HI , where the sampling process is controlled 
by a guiding function [36| . This allows to compute the 
negative-energy tails of the ground-state energy distribu- 
tion up to 12 standard deviations away from its mean. 
Similar to Ref. [37[ we consider the negative tails of the 
distribution for the reduced energy e < 1 only. We fur- 
ther use a guiding function 



G(e) = exp(a - b\e + c|") 



(2) 



in order to estimate the parameters that characterize best 
the simple sampling distributions at the three values p = 
0.83, p — pc, p = 0.86. The respective estimates are 
listed in Tab. [TTI and the guiding function for p = 0.83 is 
indicated as dashed line in Fig. [8j 

Let P(e) describe the true probability density function 
of observing a minimum-energy path with reduced energy 
e for the model under consideration (bear in mind that it 
holds that P(e) = <t e Pi(E), where e = (E- (E))/a E ). In 
order to arrive at an ImSa estimate P IS (e) that approxi- 
mates P(e), we divide the generation of the I iterations 
into two parts, consisting of J — AI and AI iterations, 
where AI is small, we consider AI = 2, 3 or 4. The gen- 
eration of the first / — AI iterations, is performed in the 
usual simple-sampling way, leading to a large (N = 10 6 ) 
pool Vi-ai- The final AI iterations should be done in a 
way that also the interesting tails of the distribution are 
sampled. For this purpose, we generate a Markov chain 
7>o S -)• V\ s -> VI s -> •■• of pools, which all are subsets 
of Vi-ai, but the sampling is done in a way that also 
the tails of the distribution of reduced energy values are 
sampled. 

The initial pool Vq S is created by picking a (uniformly 
sampled) random subset of 5 A/ edges from Vi-ai- For 
this pool now the final AI levels of the hierarchy are 
performed within one graph: The 5 A/ edges comprising 
the sampling pool can be arranged into one particular 
realization of a hierarchical graph Gai- Upon stepwise 
decimation Gai ~> Go this yields one particular edge 
with an edge energy Eq and a corresponding reduced en- 
ergy eo = (-Eq S — (E))/<te, wherein (E) and <je describe 



TABLE II: Parameters for the guiding function used during 
the importance sampling procedure for the listed values p of 
the disorder parameter. The values are obtained by fitting the 
function G(e) (see text) to the negative tail of the distribution 
of the path-energies at iteration level /. 



p 


I 


a 


b 


c 


V 


X 2 /dof 


0.83 


20 


-0.98(4) 


0.72(5) 


0.14(8) 


1.55(3) 


1.01 


Pc 


15 


-1.00(3) 


0.81(6) 


0.24(7) 


1.49(4) 


0.81 


0.86 


20 


-0.94(2) 


0.56(4) 


0.03(7) 


1.84(4) 


0.78 



the simple sampling estimate of Pj(E). Note that the 
probability in the tails is very small, hence ImSa will not 
change (E) and <je considerably. 

The Markov chain Monte Carlo step reads as follows: 
From a given sampling pool V\ we construct the next 
sampling pool V}+i using the following 2-step procedure: 

(1) randomly choose a fraction p of edges contained in 
the sampling pool Vj s and replace those edges by 
new edges chosen from the large pool Vi-ai- This 
then specifies a candidate V' for the next sampling 
pool, characterized by the reduced energy e'. 

71s 



(2) set VS.! = V with probability 



accept 



lG(e')-' 



(3) 



Set Vj^ = Vj s otherwise. 



To complete the importance sampling simulation, the 
evolution of the initial sampling pool V^ s is followed a 
number of M steps. The resulting M + 1 reduced energy 
values €q . . .€m comprise an auxiliary distribution P IS (e) 
that describes the probability by means of which a re- 
duced path energy e is visited within the ImSa procedure. 
Since the importance sampling is designed such that a 
sampling pool having reduced energy e is encountered 
with probability oc 1/G(e), the distribution P IS (e) is fur- 
ther given by the ratio P IS (e) = P(e)/G(e) (where P(e) 
describes the true pdf of observing a minimum-energy 
path with reduced energy e for the considered model sys- 
tem). As long as the guiding function G{e) provides a 
reasonable approximation to the true distribution of path 
energies, the auxiliary distribution P IS (e) obtained using 
the IS procedure is rather "flat". Thus, regarding the 
target distribution P(e) one might hope to improve on 
the negative-tail statistics provided by a simple sampling 
approach. As discussed in Ref. [Hj, successive configu- 
rations (i.e. sampling pools) encountered during an ImSa 
simulation are not independent. As a remedy one might 
consider the autocorrelation function 



X(A<) 



i+Ai) 



(Ei)(Ei 



+Ai) 



(4) 



associated to the sequence of energy values E^ obtained 
from the importance sampling procedure. The number 
of Monte Carlo steps that have to elapse until the au- 
tocorrelation function decays to 1/e gives the respective 
autocorrelation time te- Sampling pools that are sep- 
arated by r; te Monte Carlo steps can be considered 
effectively uncorrelated. Finally, the truncated sequence 
of effectively uncorrelated energy values can be analyzed 
similar to the simple-sampling data. 

The results discussed below were obtained for the 
choice p = 0.2, where we restricted the IS procedure 
to —12 < e < —1. We further performed a number of 
AI = 10 7 Monte Carlo steps to estimate the distribution 
of path energies in a target pool corresponding to / = 20. 
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FIG. 9: The main plot shows the autocorrelation functions 
x(Ai) for importance sampling simulations at p — 0.83 and 
AI — 2, 3, 4. The solid lines indicate a best fit of the re- 
spective data to a stretched exponential function x(A«) oc 
cxp{ — (Ai/rs)' 3 }. The resulting fit-parameters are listed in 
the text. The inset indicates the effective autocorrelation 
times that result from an analysis of slices of AM succes- 
sive energy values obtained during the ImSa simulation at 
AI — 3. For AM < 10 5 the autocorrelation time appears to 
increase oc AM 0,87 ' 3 ' (the corresponding fit is indicated as a 
dashed line). 



In order to assess the autocorrelation time we consid- 
ered IS simulations at a disorder parameter p = 0.83 and 
for AI = 2, 3 and 4. As it appears, the autocorrelation 
function can be well described by a stretched exponen- 
tial, i.e. x(A0 oc exp{ — (Ai/rg)^}, where j3 < 1. A 
best fit of that relation to the data yields the parame- 
ters t e » 18677,4023,1817 and (3 « 0.41,0.61,0.72 for 
AI = 2,3,4, respectively. Such a stretched exponential 
might result from a continuous sum of pure exponential 
decays [42]. Further, a stretched exponential decay of the 
energy autocorrelation function was, above the respec- 
tive critical temperature, also observed for the 2D and 
3D fully frustrated Ising model 43]. To check whether 
there are different autocorrelation times relevant on dif- 
ferent timescales of the ImSa simulation we performed 
the following analysis: we subdivided the "long" simu- 
lation run at AI = 3 into m = M/ AM "shorter" runs, 
each of length AM. The m runs of length AM are then 
considered as being independent and the characteristic 
autocorrelation time te(AM) for the shorter sequences 
is estimated. In this regard, for AM < 5 x 10 3 it was 
sufficient to consider a pure power law fit-function. For 
values of AM larger than that, a stretched exponential 
turned out to be more adequate. As shown in the inset of 
Fig.lHl we found that for AM < 10 5 the effective autocor- 
relation time is well described by an algebraic dependence 
t e {AM) = 0.24(8) x AM ' 87 ' 31 , whereas for AM > 10 5 
the value of te did not increase further, i.e. we observed 
te(AM > 10 5 ) « 4000. As pointed out above, sampling 
pools that are separated by more than te ~ 4000 Monte 



Carlo steps are effectively uncorrelated. One might now 
argue that, in order to use only uncorrelated values of 
E to construct the distribution P IS (e) (and hence P(e)), 
one should keep only every r^-th energy value. How- 
ever, in Ref. [33] the authors concluded that if during the 
M Monte Carlo steps the interval [e m i n , e m ax] is crossed 
sufficiently often, it is not necessary to discard any en- 
ergy values obtained during the IS simulation. Here, for 
AI = 3, considering the interval [—12, —1] and perform- 
ing M = 10 7 Monte Carlo steps at p = 0.83 we found a 
number of n croS s = 366 interval crossings. In agreement 
with Ref. [37J we observed that it makes no difference 
whether the sequence of energy values collected during 
the IS procedure was truncated or not, the resulting pdf 
P(e) remained almost unchanged (apart from effects due 
to the different sample sizes used to construct the pdfs). 

Below we present the results obtained for ImSa simula- 
tions at p = 0.83, p c , 0.86 considering AI — 3 and a tar- 
get distribution at / — 20 (only the simulation at p c was 
carried out for a target pool at I = 15). Once the distri- 
bution P (e) is obtained from the ImSa procedure, it can 
immediately be transformed to the desired pdf P(e). A 
comparison of the SiSa and ImSa pdf shows that the ab- 
solute probabilities in the overlapping region e S [—6, —1] 
do not coincide. This is due to the restriction of the rel- 
ative energies to the interval e G [—12,-1] during the 
ImSa simulation. One can easily account for this dis- 
crepancy by requiring that the SiSa and ImSa estimates 
coincide in the overlapping region, and by rescaling the 
ImSa estimate accordingly. In Figs. [TUT a-c). the resulting 
pdfs of observing a minimum-energy path with reduced 
energy e for p — 0.83, p c , 0.86 are shown, respectively. 
Note that the pdfs are represented using histograms that 
consist of 64 bins, each. In either case, the ImSa (SiSa) 
estimate is depicted for relative energies e < —2 (> —2). 
As evident from the figures, using the ImSa procedure 
probabilities as small as P(e) oc 10 -20 can be reached, in 
contrast to oc 10~ 6 for a SiSa approach, cf. Fig. [Sj The 
insets to Figs. [TUTa-c) indicate the normalized deviation 
of a best fit of the function G(e) — exp{a — b\e — c\ v } to 
the negative tail of P(e). Once the fit is performed the 
deviation is obtained as Agt(e) = (P(e) — G(e))/AP(e), 
where AP(e) indicates the measurement error on P(e) as 
obtained by bootstrap resampling [40| . The observation 
that Afit (considering a fit of G(e) to the ImSa data) is 



TABLE III: Exponent rj, describing the scaling of the negative 
tail of the pdf P(e), obtained by fitting the function G(e) (see 
text) to the data resulting from the ImSa procedure. From 
left to right: value p of the disorder parameter, interval over 
which the fit was performed, exponent rj (as well as 1 — 1/rj) 
and reduced chi-square \ 2 /dof . 



p 


[e_ , e^ 


-] 


V 


1-1/7) 


X 2 /dof 


0.83 


[-U, 


-1.5] 


1.42(3) 


0.30(1) 


0.88 


Pc 


[-10, 


-1] 


1.42(2) 


0.30(1) 


1.21 


0.86 


[-12, 


-1] 


1.65(6) 


0.39(2) 


1.10 
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FIG. 10: Pdf of observing a minimum- energy path with 
reduced energy e, obtained by an importance sampling Monte 
Carlo simulation in the disorder. The main plots show a semi- 
logarithmic plot of the distribution P(e) at (a) p = 0.83, (b) 
p « p c , and (c) p = 0.86. Data points at e < —2 (> —2) refer 
to the ImSa (SiSa) estimate of the pdf. In either case the error 
bars are of the order of the symbol size. The insets show the 
normalized deviation Afit of a best fit for G(e) = exp{a — 
b\x + c\ v } to the negative tail. The respective parameters rj 
are listed in Tab. Mil 



of order one and changes sign in an irregular fashion in- 
dicates that there are no systematic deviations and that 
the fit-function G(e) represents a proper approximation 
to the negative tail of the observed pdf. In particular, 
Fig. [Tola) highlights that a fit to the SiSa pdf might be 
misleading if one is interested in the true scaling behav- 
ior of P(e) as e — > — oo. For that purpose, the above 
fit-function with fitting parameters obtained for the SiSa 
data (dashed curve in the main plot) was used to com- 
pute the normalized deviation to the ImSa data (dashed 
line in the inset). Referring to this one finds rather strong 
systematic deviations where, e.g., Afi t (— 11) w 10. The 
parameters 77 that correspond to a best fit to the ImSa 
data are listed in Tab. Mil 

As pointed out above, previous studies suggested that 
the exponent rj is related to the energy fluctuation ex- 
ponent uj (listed in Tab. UJ by means of the expression 
uj = 1 — 1/77. Here, for the minimum- energy path problem 
on hierarchical lattice graphs we find that this expression 
holds for all values of p thus considered. To facilitate 
comparison, the values 1 — 1 / 77 are listed in Tab. Hill 



IV. CONCLUSIONS 

In the presented article we have investigated a particu- 
lar MWP problem on MK hierarchi cal g raphs. It is quite 
similar to earlier polymer 12, H il| and optimal 

path problems [fl l2ll l32l |33|. The important difference 
is that for a considerable fraction of negative edge ener- 
gies, the total energy of a path may be reduced by tak- 
ing longer paths (which leads to a different universality 
class). In the same fashion as the optimal path problem 
on hierarchical graphs, studied in Ref. [2l|, corresponds 
to the (generic) non-directed optimal path problem [loj . 
the minimum-energy path problem studied here corre- 
sponds to the negative-energy percolation problem [l3[ 
in which there is a path forced onto the system and where 
the disorder is "weak" enough to render the appearance 
of loops irrelevant for the scaling behavior of the path. 

Here, the scaling properties of the MWP obtained af- 
ter the decimation of huge hierarchical graphs change 
with increasing edge-disorder, leading from a phase where 
the path displays a self-affine scaling behavior to a phase 
where the path displays a (statistically) self-similar scal- 
ing behavior. We characterized the respective phase tran- 
sition by monitoring the length and energy of the MWPs 
as function of the disorder and quantified the scaling 
behavior of the observables (and their fluctuations) by 
means of proper critical exponents. While the scaling 
of the observables off criticality can be explained intu- 
itively, the scaling behavior found at the critical point of 
the model is nontrivial and compares well to the scaling 
observed for the optimal path problem on the same hi- 
erarchical lattice in the limit of "strong" disorder [2~lj |. 
However, note that the precise optimization criteria of 
both models are slightly different: while in the optimal 
path problem [ToL l2lj one aims to minimize the largest 
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energy along a sub-path, here one strives after minimiz- 
ing the sum of energies along a sub-path. Note that this 
was already realized for the respective models on regular 
lattice graphs, where quite similar scaling exponents for 
the average path length were found in dimensions d = 2 
through 6, see Ref. [la ]. 

Further, we performed an importance sampling simula- 
tion for the ground-state energy distribution of the paths 
and confirmed that it is consistent with a Tracy- Widom 
scaling form, similar to the directed polymer in a random 
medium [13] )• Using the importance sampling procedure 
allowed for an analysis of the ground-state energy dis- 
tribution down to probabilities as small as oc 10~ 20 (in 
contrast, a convenient simple sampling approach only al- 
lows to reach oc 10~ 6 ), i.e. up to 12 standard deviations 
away from the mean of the respective distribution. This 
leads to a precise estimate of the scaling behavior of the 



negative tail of the ground-state energy distribution. For 
all values of the disorder parameter considered here, the 
respective exponents rj could be related to the energy- 
fluctuation exponents lo via the relation lu = 1 — l/i]. 



Acknowledgments 

OM acknowledges financial support from the DFG 
(Deutsche Forschungsgemeinschaft) under grant 
HA3 169/3-1. The simulations were performed at 
the HPC Cluster HERO, located at the University of 
Oldenburg (Germany) and funded by the DFG through 
its Major Instrumentation Programme (INST 184/108-1 
FUGG) and the Ministry of Science and Culture (MWK) 
of the Lower Saxony State. 



K. Kremer, Z. Phys. B 45, 149 (1981). 
M. Kardar and Y. C. Zhang, Phys. Rev. Lett. 58, 2087 [22 
(1987). 

B. Derrida, Physica A 163, 71 (1990). [23 
P. Grassberger, J. Phys. A 26, 1023 (1993). 
R. Parshani, L. A. Braunstein, and S. Havlin, Phys. Rev. [24 
E 79, 050102 (2009). [25 
M. Cieplak, A. Maritan, and J. R. Banavar, Phys. Rev. [26 
Lett. 72, 2320 (1994). 

O. Melchert and A. K. Hartmann, Phys. Rev. B 76, [27 
174411 (2007). 

O. Melchert and A. K. Hartmann, Comp. Phys. Comm. [28 
182, 1828 (2011). [29 
K. Schwarz, A. Karrenbauer, G. Schehr, and H. Rieger, [30 
J. Stat. Mech. 2009, P08022 (2009). 

N. Schwartz, A. L. Nazaryev, and S. Havlin, Phys. Rev. [31 
E 58, 7642 (1998). 

T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, [32 
Introduction to Algorithms, 2nd edition (MIT Press, 
2001). [33 
O. Melchert and A. K. Hartmann, Phys. Rev. B 79, 
184402 (2009). [34 
O. Melchert and A. K. Hartmann, New. J. Phys. 10, [35 
043039 (2008). 

L. Apolo, O. Melchert, and A. K. Hartmann, Phys. Rev. [36 
E 79, 031103 (2009). 
O. Melchert, L. Apolo, and A. K. Hartmann, Phys. Rev. [37 
E 81, 051108 (2010). 

O. Melchert, A. K. Hartmann, and M. Mezard, Phys. [38 
Rev. E 84, 041106 (2011). 

C. Norrenbrock, O. Melchert, and A. K. Hartmann [39 
(2012), preprint: arXiv:1205.1412. 

G. Claussen, L. Apolo, O. Melchert, and A. K. Hart- 
mann, Phys. Rev. E 86, 056708 (2012). [40 
R. K. Ahuja, T. L. Magnanti, and J. B. Orlin, Network 
Flows: Theory, Algorithms, and Applications (Prentice [41 
Hall, 1993). 

D. Stauffer and A. Aharony, Introduction to Percolation [42 
Theory (Taylor and Francis, London, 1994). [43 
M. Cieplak, A. Maritan, M. R. Swift, A. Ahattacharya, 
A. L. Stella, and J. R. Banavar, J. Phys. A 28, 5693 



(1995). 

O. R. Salmon, B. T. Agostini, and F. D. Nobre, Phys. 
Lett. A 374, 1631 (2010). 

R. Teodoro, C. G. Bezerra, A. M. Mariz, and F. A. 
da Costa, Preprint: arXiv:1005.3863vl (2010). 

A. J. Bray and S. Feng, Phys. Rev. B 36, 8456 (1987). 
S. Boettcher, Eur. Phys. J. B 33, 439 (2003). 

B. Derrida and R. B. Griffiths, Europhys. Lett. 8, 111 
(1989). 

P. Doussal and J. Machta, J. Stat. Phys. 64, 541578 
(1991). 

P. Devillard, Phys. Rev. Lett. 70, 1124 (1993). 

M. S. Cao, J. Stat. Phys. 71, 51 (1993). 

Y. Shussman and A. Aharony, J. Stat. Phys. 80, 147 

(1995). 

C. Monthus and T. Garel, Phys. Rev. E 77, 021132 
(2008). 

A. Hansen and J. Kertesz, Phys. Rev. Lett. 93, 040601 
(2004). 

S. V. Buldyrev, S. Havlin, and H. E. Stanley, Phys. Rev. 
E 73 (2006). 

A. K. Hartmann, Phys. Rev. E 65, 056102 (2002). 

A. Engel, R. Monasson, and A. K. Hartmann, J. Stat. 

Phys. 117, 387 (2004). 

M. Korner, H. G. Katzgraber, and A. K. Hartmann, J. 
Stat. Mech. p. P04005 (2006). 

C. Monthus and T. Garel, Phys. Rev. E 74, 051109 
(2006). 

S. Wolfsheimer, B. Burghardt, and A. K. Hartmann, Al- 
gorithms for Molecular Biology 2, 9 (2007). 
W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. 
Flannery, Numerical Recipes in C (Cambridge University 
Press, Cambridge, U.K., 1992). 

A. K. Hartmann, Practical Guide to Computer Simula- 
tions (World Scientific, Singapore, 2009). 
T. Halpin-Healy and Y.-C. Zhang, Phys. Rep. 254, 215 
(1995). 

D. C. Johnston, Phys. Rev. B 74, 184430 (2006). 

G. Franzese, A. Fierro, A. De Candia, and A. Coniglio, 
Physica A 257, 376 (1998). 



